Advanced numerical engineering

& Gaussian Process

Y. Richet (yann.richet@irsn.fr)

credits: N. Durrande, R. Le Riche, V. Picheny, N. Garland

2018

…

General metamodel engineering approach
Optimization basics & engineering context
Inversion basics & engineering context
Improved numerical engineering algorithms & models

Overview

What was/is engineering ?

Before 70's
Analytical approach to roughly design nuclear plants, rockets, cars, …

The Unreasonable Effectiveness of Mathematics, E. Wigner

Since 70-80's Experimental/numerical engineering to simulate (vs. solve) fine multi-physics, unstable systems, …

"I am an old man now, and when I die and go to heaven there are two matters on which I hope for enlightenment.
One is quantum electrodynamics, and the other is the turbulent motion of fluids.

And about the former I am rather optimistic.", H. Lamb

Engineering: industry & standards

Engineering: industry vs. standards

Industrial products:

wide product target identification
optimization of products profitability
optimization of investments

Regulatory control:

define industrial-field standards
- end-user safety
- ecological & resources parcimony
- (fair market)
check products conformity
evolve standards to best knowledge

Optimization (reminder)

f^{*} = m i n_{x \in S \subset R^{d}} f (x)

$f^* = min_{x\ \in\ S\ \subset\ \mathbb{R}^d} \ \ f(x)$

x^{*} = a r g m i n_{x \in S \subset R^{d}} f (x)

$x^* = argmin_{x\ \in\ S\ \subset\ \mathbb{R}^d} \ \ f(x)$

$S$ is the input domain for $x$
$d$ is the dimension of $x$
$f$ is the cost function (ex. stress, risk, power, temperature, …)

for sake of simplicity and by convention, we consider minimization

Optimization

Optimization - basics

Gradient / Newton methods

Optimization - basics

Evolutionnary algorithms (CMA-ES, PSO, …)

Optimization - basics vs. engineering

Optimization in engineering context raises many issues:

one $f$ evaluation is expensive (CPU/time/user)
dimension of problem input $d$ may be large
time to result is constrained
added-value from optimization is initially unknown

But basic methods are not suitable:

Gradient / Newton methods
- gradients have to be computed (finite differences ?)
- local optimization, while we need global
- handle noise of $f$ evaluation ?
Evolutionnary algorithms
- too much calls of $f$

Optimization - basics vs. engineering

Optimization in engineering context raises many issues:

one $f$ evaluation is expensive (CPU/time/user)
dimension of problem input $d$ may be large
time to result is constrained
added-value from optimization is initially unknown

But basic methods are not suitable

So, true practice is often very rough:

change one parameter at a time (to avoid mistakes)
$f$ evaluations performed should be enough different
startup with some assumptions on the solution
- intrinsinc physics
- consider monotonous effects (even approximate)

Optimization - engineering

One-at-a-time optimization

Optimization - engineering++

Bayesian optimization (EGO,…)

Inversion

{x^{*}} = a r g_{x \in S \subset R^{d}} {f (x) = T}

$\{ x^* \} = arg_{x\ \in\ S\ \subset\ \mathbb{R}^d} \ \{ f(x) = T \}$

{x^{*}} = a r g_{x \in S \subset R^{d}} {f (x) < T}

$\{ x^* \} = arg_{x\ \in\ S\ \subset\ \mathbb{R}^d} \ \{ f(x) < T \}$

$S$ is the input domain for $x$
$d$ is the dimension of $x$
$f$ is the cost function (ex. stress, risk, power, temperature, …)

Inversion

Numerical engineering algorithms…

uncertainties propagation:
- "What spreading of output when inputs follow given random laws ?"
- random sampling (~ Monte Carlo), dedicated models (polynomial chaos)
sensitivity analysis
- "What share of output spreading is due to random inputs ?"
- relative indices: $V[S[Y|X]]/...$ (Sobol, HSIC, …), random sampling, dedicated models (Fourier, polynomial chaos)
parameters screening
- "What output behaviour for each input ?"
- screening designs (sparse), OaT designs (Morris)
optimization
- "What worst/best value possible for output ?"
- bayesian optimization (EGO)
inversion
- "What input values give some output value ?"
- bayesian inversion (Bichon's, TIMSE, SUR)

Metamodel based engineering

Basic idea

f

$\color{red} f$

Basic idea

To limit number of $\color{red} f$ evaluations:

create a $\color{blue}{model\ of\ f}$ based on few evaluations $x={X}$

Basic idea

To limit number of $\color{red} f$ evaluations:

create a $\color{blue}{model\ of\ f}$ based on few evaluations $x={X}$

But it is risky to take decision only based on a single model …

… However, we would expect the minimum to be not so far from model's one.

Basic idea

To limit number of $\color{red} f$ evaluations:

create a $\color{blue}{model\ of\ f}$ based on few evaluations $x={X}$

But it is risky to take decision only based on a model …

=> Diversifiy the model

… However, we would expect the minimum to be not so far from model's one.

=> Take model minimum just as a clue

Basic idea

To limit number of $\color{red} f$ evaluations:

create a $\color{blue}{model\ of\ f}$ based on few evaluations $x={X}$

But it is risky to take decision only based on a model …

=> Diversifiy the model: [EXPLORE]

… However, we would expect the minimum to be not so far from model's one.

=> Take model minimum just as a clue: [EXPLOIT]

Basic idea

To limit number of $\color{red} f$ evaluations:

create many $\color{blue}{models\ of\ f}$ based on few evaluations $x={X}$

Conditional Gaussian Process

Suitable to "interpolate" few evaluations
Suitable to sample in a model family

Conditional Gaussian Process

Suitable to "interpolate" few evaluations
Suitable to sample in a model family
Give an explicit (Gaussian) density marginal on x

Conditional Gaussian Process

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

Model pointwise mining

define a suitable criterion for engineering target
find its maximum over $X$
sample $f$ at this new point

Optimization

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

To trade-off between exploration/exploitation we will consider:

distribution of $M(x)$ : [EXPLORE]
set of $x$ where $M(x) < min \{ f(X) \}$ : [EXPLOIT]

Efficient Global Optimization

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

Let's define the Probability of Improvement:

P I (x) = P [M (x) < m i n {f (X)}]

$\color{violet}{PI}(x) = P[M(x) < min \{ f(X) \}]$

which is analytical thanks to $M$ properties…

EGO (reminder)

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

u (x) = \frac{m i n {f (X)} - m (x)}{s (x)}

$u(x) = {{min\{f(X)\} - m(x)} \over {s(x)}}$

P I (x) = p_{N} (u (x))

$\color{violet}{PI}(x) = p_{\mathcal{N}} \big( u(x) \big)$

EGO (reminder)

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

P I (x) = P [M (x) < m i n {f (X)}]

$\color{violet}{PI}(x) = P[M(x) < min \{ f(X) \}]$

But, $x$ for highest $PI$ is often close to best $X$ …

EGO (reminder)

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

Let's define the Expected Improvement:

E I (x) = E [(m i n {f (X)} - M (x))^{+}]

$\color{violet}{EI}(x) = E[(min \{ f(X) \} - M(x))^+]$

which is (also) analytical thanks to $M$ properties…

EGO (reminder)

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

u (x) = \frac{m i n {f (X)} - m (x)}{s (x)}

$u(x) = {{min\{f(X)\} - m(x)} \over {s(x)}}$

E I (x) = s (x) (u (x) p_{N} (u (x)) + d_{N} (u (x)))

$\color{violet}{EI}(x) = s(x)\ \big(\ u(x)p_{\mathcal{N}}(u(x)) + d_{\mathcal{N}}(u(x))\ \big)$

EGO (reminder)

"Efficient Global Optimization of Expensive Black-Box Functions"- Jones, Schonlau, Welch,
(Journal of Global Optimization, December 1998)

EGO:
Maximize $EI(x)$ (*), compute $f$ there, add to $X$ .
Repeat until …

+ good trade-off between exploration and exploitation
+ requires few evaluations of $f$
- often lead to add close points to each others …
Which is not very comfortable for kriging numerical stability
- "one step lookahead" (myopic) strategy
- rely on model suitability to $f$

(*) using standard optimization algorithm: BFGS, PSO, DiRect, …

EGO (reminder)

Support for parallel evaluations of $f$

Obviously major feature of optimization:
- make profit of HPC investment
- even multi-CPU standalone computers
- non-parallel numerical solvers become relevant (again)

EGO++:

GP-consistency $q-EI$ criterion (quite hard)
Batch-sequential heuristic "Constant Liar":
1. instead of evaluating $f(x^*)$ , assume $f(x^*) = min\{f(X)\}$
2. update (including $(x^*,min\{f(X)\})$ in $X$ ) and optimize EI
3. …

Inversion

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

{x^{*}} = a r g_{x \in S \subset R^{d}} {f (x) < T}

$\{ x^* \} = arg_{x\ \in\ S\ \subset\ \mathbb{R}^d} \ \{ f(x) < T \}$

Build your own inversion criterion …

what [approx] trend(x) will be useful ?
draw $s(x)$ , $m(x)$ , $T-m(x)$ , …

Build your own inversion criterion …

Efficient Global Inversion - Bichon's

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

Let's define the Expected Excursion:

E E (x) = E [(s_{n} (x) - | T - M (x) |)^{+}]

$\color{violet}{EE}(x) = E[(s_n(x) - |T - M(x)|)^+]$

which is (also) analytical thanks to $M$ properties…

EGI - Bichon's

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

v (x) = \frac{T - m (x)}{s (x)}

$v(x) = {{T - m(x)} \over {s(x)}}$

\begin{aligned} E E (x) = s (x) ( & p_{N} (v (x) + 1) - p_{N} (v (x) - 1) \\ - v (x) (2 p_{N} (v (x)) - p_{N} (v (x) + 1) - p_{N} (v (x) - 1)) \\ - (2 d_{N} (v (x)) - d_{N} (v (x) + 1) - d_{N} (v (x) - 1)) \\ ) \end{aligned}

$\small{ \begin{align*} \color{violet}{EE}(x) = s(x)\ \big(\ & \ p_{\mathcal{N}}(v(x)+1) -p_{\mathcal{N}}(v(x)-1) \\ &- v(x) (2p_{\mathcal{N}}(v(x)) - p_{\mathcal{N}}(v(x)+1) -p_{\mathcal{N}}(v(x)-1)) \\ &- ( 2d_{\mathcal{N}}(v(x)) - d_{\mathcal{N}}(v(x)+1) -d_{\mathcal{N}}(v(x)-1) ) \ \\ &\big) \end{align*} }$

EGI - Bichon's

"Sequential design of computer experiments for the estimation of a probability of failure" - Bect, Ginsbourger, Li, Picheny, Vazquez, (Statistics and Computing 2012)

EGI:
Maximize $EE(x)$ (*), compute $f$ there, add to $X$ .
Repeat until …

+ trade-off between exploration and exploitation
+ requires few evaluations of $f$
- often lead to add close points to each others …
Which is not very comfortable for kriging numerical stability
- "one step lookahead" (myopic) strategy
- rely on model suitability to $f$

(*) using standard optimization algorithm: BFGS, PSO, DiRect, …

EGI - Bichon's

Model conditional uncertainty mining

define a suitable criterion for engineering target
compute its integral conditionally to any $x$ added
find the minimizer of $x$
sample $f$ at this new point

Inversion

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

Exploration >> exploitation, which is interesting for identification problems, where discrete solution never apply (like inversion).

Stepwise Uncertainty Reduction

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

Let's define the Probability of Excursion:

P E (x) = P [M (x) < T]

$\color{violet}{PE}(x) = P[M(x) < T]$

… and the Uncertainty of Excursion:

S E (x) = P [M (x) < T] \times P [M (x) > T]

$\color{violet}{SE}(x) = P[M(x) < T] \times P[M(x) > T]$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

v (x) = \frac{T - m (x)}{s (x)}

$v(x) = {{T - m(x)} \over {s(x)}}$

S E (x) = p_{N} (v (x)) \times (1 - p_{N} (v (x)))

$\color{violet}{SE}(x) = p_{\mathcal{N}} \big( v(x) \big) \times (1-p_{\mathcal{N}} \big( v(x) \big))$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

y_{n e w} \sim N (m (x_{n e w}), s^{2} (x_{n e w}))

$y_{new} \sim \mathcal{N}(m(x_{new}),s^2(x_{new}))$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

y_{n e w} \sim N (m (x_{n e w}), s^{2} (x_{n e w}))

$y_{new} \sim \mathcal{N}(m(x_{new}),s^2(x_{new}))$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

y_{n e w} \sim N (m (x_{n e w}), s^{2} (x_{n e w}))

$y_{new} \sim \mathcal{N}(m(x_{new}),s^2(x_{new}))$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

y_{n e w} \sim N (m (x_{n e w}), s^{2} (x_{n e w}))

$y_{new} \sim \mathcal{N}(m(x_{new}),s^2(x_{new}))$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

y_{n e w} \sim N (m (x_{n e w}), s^{2} (x_{n e w}))

$y_{new} \sim \mathcal{N}(m(x_{new}),s^2(x_{new}))$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

y_{n e w} \sim N (m (x_{n e w}), s^{2} (x_{n e w}))

$y_{new} \sim \mathcal{N}(m(x_{new}),s^2(x_{new}))$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

y_{n e w} \sim N (m (x_{n e w}), s^{2} (x_{n e w}))

$y_{new} \sim \mathcal{N}(m(x_{new}),s^2(x_{new}))$

SUR

M (x) = N (m (x), s^{2} (x))

$\color{blue}M(x) = \mathcal{N}(m(x),s^2(x))$

$m(x) = C(x)^T C(X)^{-1} f(X)$
$s^2(x) = c(x) - C(x)^T C(X)^{-1} C(x)$
$C$ is the covariance kernel $C(.) = C(X,.)$ , $c(.) = C(x,.)$

y_{n e w} \sim N (m (x_{n e w}), s^{2} (x_{n e w}))

$y_{new} \sim \mathcal{N}(m(x_{new}),s^2(x_{new}))$

Reference

KrigInv: An efficient and user-friendly implementation of batch-sequential inversion strategies based on Kriging

Clément Chevalier, IMSV University of Bern
Victor Picheny, INRA Toulouse
David Ginsbourger, IMSV University of Bern

Implements

Pointwise criterions: Bichon, Ranjan, TMSE
SUR criterions: TIMSE, SUR, Jn

Conclusions

Pointwise vs. conditional uncertainty

Instead of a local criterion (like $SE(x)$ ),
use a global gain which integrates local criterion:

S E (x) = P [T < M_{n} (x)] \times P [T > M_{n} (x)] ↓ S U R (x) = \int_{S} P [T < M_{n + {x}} (y)] \times P [T > M_{n + {x}} (y)] d y

$SE(x) = P[T < M_n(x)] \times P[T > M_n(x)] \\ \downarrow \\ SUR(x) = \int_S P[ T < M_{n+\{x\}}(y)] \times P[ T > M_{n+\{x\}}(y)] \ dy$

i.e. We search for $x$ which, once added to $X$ , most reduces $\int_S SE$ .

Pointwise vs. conditional uncertainty

Instead of a local criterion (like $EI(x)$ ),
use a global gain which integrates local criterion:

E I (x) = E [m i n {f (X_{n})} - M_{n} (x)]^{+} ↓ I E C I (x) = \int_{S} E [m i n {f (X_{n} \oplus x)} - M_{n + {x}} (y)]^{+} d y

$EI(x) = E[min \{ f(X_n) \} - M_n(x)]^+ \\ \downarrow \\ IECI(x) = \int_S E[min \{ f(X_{n} \oplus x) \} - M_{n+\{x\}}(y)]^+\ dy$

i.e. We search for $x$ which, once added to $X$ , most reduces $\int_S EI$ .

Stepwise Uncertainty Reduction

SUR criterions are:

more flexible than local ones () to handle:
- constraints:
  $m i n_{x \in S \subset R^{d}} {f (x) : g (x) \leq 0}$ $min_{x\ \in\ S\ \subset\ \mathbb{R}^d} \ \{ f(x) : g(x) \le 0\}$
- robustness:
  $m i n_{x \in S \subset R^{d}} {f (x) : g (x, u) \leq 0, \forall u \in U}$ $min_{x\ \in\ S\ \subset\ \mathbb{R}^d} \ \{ f(x) : g(x,u) \le 0, \forall u \in U \}$
but much more costly to evaluate ( $\int_S$ )

Constrained optimization

Basically using $EI$ :

E I (x) = 1_{g (x) \leq 0} E [m i n {f (X)} - M (x)]^{+}

$EI(x) = \mathbf{1}_{g(x) \le 0} \ E[min \{ f(X) \} - M(x)]^+$

… but this excludes to sample $x$ if $g(x)>0$ .

Or using $IECI$ :

I E C I (x) = \int_{S} 1_{g (y) \leq 0} E [m i n {f (X_{n} \oplus x)} - M_{n + {x}} (y)]^{+} d y

$IECI(x) = \int_S \mathbf{1}_{g(y) \le 0} \ E[min \{ f(X_{n} \oplus x) \} - M_{n+\{x\}}(y)]^+ \ dy$

Tips & tricks

On-shelf available (tech.)

General optimization

Languages:

R, Python, Matlab/Octave
C++, Java, C, Fortran

Libraries:

NLOpt (C, C++, Fortran, Matlab/Octave, Python, Julia, R)
DiRect, CRS, COByLA, Nelder-Mead, BFGS, …
SciPy.optimize (Python) BFGS, COByLA, Nelder-Mead, …
optim (R) Nelder-Mead, BFGS, Simulated-Annealing, …

State-of-the-Art algorithms: EGO, CMA-ES, NSGA2, PSO
(C, C++, Fortran, Matlab/Octave, Python, R, Scilab, Julia, …)

On-shelf available (tech.)

Kriging

Languages:

R, Python, Matlab/Octave
C++, Java,

Libraries:

DiceKriging/KerGP (R)
SKL (Matlab/Octave)
geoR (R)
GPy/GPflow (Python)
OpenTURNS (Python)

If problem remains too hard…

Go deeper in $f$ .

Get more runs
- numerical model: more CPU, HPC, cloud
Get more prior informations
- boundary conditions
- trending hypothesis
- simplified model (multi-level opt.)
Get derivatives
- physical model: intrinsic physics
- numerical model: integrated|auto differentiation
  - source code differentiation
  - dedicated framework (ex. TensorFlow)